June 2003 



Successive Umbrella Sampling 

Peter VirnarQ and Marcus Miiller 
Institut fiir Physik, WA331, Universitat Mainz, Staudinger Weg 7, 55099 Mainz, Germany 

(Dated: February 2, 2008) 

We propose an extension of umbrella sampling in which the pertinent range of states is subdivided 
in windows that are sampled consecutively and linked together. Extrapolating results from one 
window we estimate a weight function for the neighboring window. We present a detailed analysis 
and demonstrate that the error is controlled and independent from the window sizes. The analysis 
also allows us to detect sampling difficulties. The efficiency of the algorithm is comparable to a 
multicanonical simulation with an ideal weight function. We exemplify the computational scheme 
by simulating the liquid- vapor coexistence in a Lennard-Jones system. 
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Calculating and overcoming free energy barriers re- 
mains a considerable challenge in computational con- 
densed and soft matter physics. Applications are many- 
fold and range from protein folding |lj over quantum sys- 
tems |2( to the study of first order phase transitions [3J, 
nucleation or glassy systems |4|. 

Various sophisticated sampling techniques have been 
devised to estimate free energy differences. Though 
methods are general we shall use the language of a liquid- 
vapor transition, where the number of particles n is 
the order parameter of the transition. In this case the 
Helmholtz free energy A can be computed from the prob- 
ability distribution P[n]: A[n] = —k^,Th\P[n] + const., 
where fee stands for Boltzmann's constant and T for tem- 
perature. Free energy barriers correspond to regions of 
extremely low probability, which make efficient sampling 
difficult. 

Multicanonical methods [5j modify the Hamiltonian 
in order to sample a range of densities uniformly. To 
this end, one adds a weight function w[n] to the orig- 
inal Hamiltonian such that the simulated distribution 
fsimM = P[ n ] exp(— w[n]) becomes fiat for the choice 
of w[n] w In P[n]. Unfortunately, P[n] is a priori un- 
known, and several methods have been explored to esti- 
mate w[n\: (i) Histogram-reweighting techniques 6] alle- 
viate this problem by performing a sequence of weighted 
simulations and extrapolations starting at a point where 
barriers are small and the system explores a broad range 
of n. More sophisticated methods combine results of 
multiple histograms Q. (ii) Multicanonical Recursion 
H conducts a series of short trial runs. After each run 
w[n] is adjusted until the simulation can access all rele- 
vant states. The weight factors can also be self-adjusted 
during the simulation 0,EJE1E]]- However, detailed 
balance 0] is violated in this process and separation of 
statistical and systematic errors becomes difficult, (hi) 
Weight factors can also be obtained from the transition 
probabilities between macrostates 0, Q, E| . 

In umbrella sampling the pertinent range of 

macrostates, defined by the number of particles n, is sub- 
divided into m overlapping windows of length uj. In this 



work we propose an extension in which adjacent windows 
are sampled consecutively. Successive umbrella sampling 
offers substantial advantages: It allows us to generate 
weight factors for subsequent windows by means of ex- 
trapolation. In contrast to histogram-reweighting, simu- 
lations can be performed anywhere in the phase diagram 
without prior knowledge of a weight function. In contrast 
to self-adjusting methods 0, El El , w[n] stays fixed 
during the run and detailed balance is not violated. Al- 
though these powerful schemes are suitable for generating 
weight functions, they require in principle an additional 
multicanonical production run whenever errors need to 
be determined exactly. In our approach, every run con- 
tributes uniformly to statistics. Additionally, it does not 
involve any adjustable parameters (e.g., the modification 
factor / in Ref. ||). 

In the following sections we demonstrate the method 
by simulating the liquid-vapor coexistence in the fiVT- 
ensemble. Beads interact via a truncated and shifted 
Lennard-Jones potential with parameters e and a: 
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The box length of our reference system is 6.74 cr, tem- 
perature T — 0.85437 4. Fig. n shows the probability 
distribution as a function of particle number n at coex- 
istence. The peak at low density corresponds to the gas 
phase, the peak at high density to the liquid phase. Both 
are separated by a free energy barrier (here: 15.2 k&T). 

The guiding idea of successive umbrella sampling is 
to investigate one small window after the other, starting 
at zero density. A histogram Hk[n] monitors how often 
each state is visited in the fcth window [ku, (k + 
Let Hki = Hk[ku>] and Hk r = Hk[{k + l)w] denote the 
values of the kth. histogram at its left and right bound- 
ary, respectively, and = Hk r /Hki characterize their ra- 
tio. After a predetermined number of insertion/deletion 
Monte Carlo (MC) moves per window, the (unnormal- 
ized) probability distribution can be estimated recur- 
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FIG. 1: P[n] for a LJ system at T=0.85437 f = 0.85 T c , L = 
6.74 a. The system is at coexistence because the area below 
both peaks is equal [l7||. Else, coexistence can be established 
by: PcocxN = P[n] ■ exp((^ cocx - fi)nk B T). Dashed lines: 
fsim[tt] from a multicanonical simulation with a very good 
weight function w[n]. Inset: free energy barrier. 
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when n G [koj, (k + l)ui]. Probability ratios in Eq. @ 
correspond to free energy differences. Care has to be 
taken at the boundaries of a window to fulfill detailed 
balance flif . 

As we are sampling one window after the other, effi- 
ciency can be increased by combining the algorithm with 
the multicanonical concept. In a weighted simulation we 
replace H[n] in Eq. (j5J) by H[n] exp(u>[n]). A good esti- 
mate for the weight function may be obtained by extrap- 
olation: After P[n] is determined according to Eq. J2J, 
w[n] = ln(P[n]) is extrapolated quadratically into the 
next window. The first window is usually unweighted. 
If in this case states are not accessible, w[n] might be 
altered by a fixed number of k#T in each iteration step. 

Let us investigate how the computational effort de- 
pends on the choice of the window size. It has been 
suggested in the literature 0, 0] that small windows 
reduce computational effort by a factor of m: if r is the 
time to obtain a predetermined statistical error Ai in a 
single window, then t oc lo 2 / Ai and computation time 
i C pu oc tow 2 /Ai. In the limit of a single large window, 



m'=l, this yields t' cpn = (moj)* 
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In Fig. we 



compare the error which was obtained from multicanon- 
ical runs with different window sizes but an equal total 
number of Monte Carlo steps 01 ■ All distributions were 
normalized to unity before the error was calculated. As 
expected, errors are evenly distributed for low and in- 
termediate densities. For large densities deviations occur 
and the relative error becomes larger. The absolute value 
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FIG. 2: Comparison of (normalized) weighted umbrella sam- 
pling runs with different window sizes and an average number 
of 2 million MC steps per state : (a) one section runs with 
window size ui = 218, (b) m = 11 section runs with window 
size lu = 20, (c) w = 1 (i.e., two states per window). The 
precision of the factors exp(«j(n)), which have been used to 
generate the data, is better than 2% for (a) and (b) and 30% 
for (c). Errors correspond to single runs and were determined 
as the standard deviation over 400 runs. 



of P[n], however, is small in this region (cf. Fig.^) such 
that the error does not influence normalization. Con- 
trary to the simple argument above the error is roughly 
independent of the window size to. To simplify the discus- 
sion, one can also set P[Q] = 1 and let errors accumulate 
from low to high densities (cf. Fig.[3J). Using Eq. J2J and 
assuming that there are no correlations between neigh- 
boring intervals we obtain 
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If errors are comparable for all windows, the overall error 
for m windows is of order 0(y/m) bigger as for m = 1. 
To compensate this, sampling time in each window has 
to be increased by a factor of m. Hence, computational 
effort for a given error is always of order (raw) 2 = A^ 2 and 
does not depend on the number of windows into which 
the sampling range is subdivided. We conclude that suc- 
cessive umbrella sampling is as fast as a multicanonical 
simulation with a very good weight function. Note fur- 
ther that a rather inaccurate estimate for w[n] leads only 
to a slight increase in the error (see figure caption in 
Fig. [21). In the following we focus on the smallest win- 
dow size u) = 1. One should keep in mind, however, that 
weighting in a single variable might not always be suffi- 
cient when free energy landscapes become more complex 

Subsequently, we present a detailed error analysis of 
our method. Eq. only accounts for statistical errors 
within an interval and neglects correlations between ad- 
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FIG. 3: Comparison of accumulated statistical error (Eq. £J) 
and independently measured error for (gi=l): (a) 4000 runs 
with 0.5 million MC steps per state |lfif . (b) 400 runs with 
2 million MC steps per state. Measured and accumulated 
error cannot be distinguished. Inset: acceptance rate of in- 
sertion/deletion attempts. 



FIG. 4: Relative systematic errors for a run with 0.5 million 
MC steps per state (compare with Fig. 0). Inset: Relation 
between relative systematic and statistical error. Line cor- 
responds to Eq. ©, o denotes the relations determined by 
simulation. 



jacent windows and possible systematic errors when es- 
timating probability ratios. The former may occur, be- 
cause we use the ending configuration of the kth window 
as the starting configuration of the (k + l)th window. 
Accumulating the relative statistical errors of individual 
ratios Ar^ in Eq. 10 , and comparing the result with the 
independently measured statistical error, we can gauge 
the importance of correlations between neighboring win- 
dows. In Fig.|3two features can be identified. First, large 
errors are obtained for high densities because it becomes 
more and more difficult to insert particles (cf. Fig. [3] in- 
set). Secondly, we find good agreement between both 
curves if the number of Monte Carlo steps per state is 
large. If the computational effort is reduced, curves only 
match for low and medium densities because in these re- 
gions correlations decay fast. Although this test is not 
crucial as errors can be measured directly, it enables us to 
check if high density configurations are well equilibrated. 
If one is only interested in a rough estimate, the total er- 
ror of a simulation can be preset before sampling starts. 
For this purpose, the termination condition is changed 
from a given number of Monte Carlo steps to a prede- 
fined error Ai in each window. Then, the total error 
adds up to y/rn ■ Ai according to Eq. J2J 21] . 

In addition to the statistical error, the use of Eq. (J2J) 
imparts a small systematic error onto the probability dis- 
tribution because the average of a probability ratio is 
larger than the ratio over two probability averages. An 
expansion of the former leads to 
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(Hi r ). (H^ I H^) corresponds to the probability ratio 
as it is measured in simulation whereas (Hi r )/(Hi r ) is 
equal to P[(i+l)w]/P[iw]. It is therefore advantageous to 
determine the probability ratio as the ratio of the sums 
of H^ and Hu over all runs (and not as the sums over 
the ratios). The last factor is typically larger than unity 
because fluctuations of Hu and Hi r are anticorrelated. 

For the special case u = 1, Hu + H r = H where H 
denotes the total number of entries in the ith histogram. 
In this case, the last factor in Eq. Q takes the form 
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In the context of a grandcanonical simulation, the sys- 
tematic error corresponds to a shift of the chemical po- 
tential of the order fcsTA^ 1 ^. Comparing this value 
with the relative statistical error Ar^ of the ratio, we ob- 
tain 
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where SHu = Hu — (Hu) and similarly SHi r = Hi r 



Relation © allows us to compute the systematic er- 
ror from the known statistical error. To test the ap- 
proximation we determined an accurate weight func- 
tion for P(n=14)/P(n=13) (gas peak) by a long run. 
(H (li) / H (13)) was computed by averaging over a large 
number of short runs with #=30, 50, 80, 200, 800 and 
2400 MC steps, respectively. The true systematic offset 
is the deviation of this ratio from 1. From the results 
in the inset of Fig. 0] we conclude that Eqs. (00) yield 
an excellent approximation for the systematic error, even 
when errors become large. For u> > 1 we expect that the 
systematic error of a single ratio is also of the order of 

Ar 2 . 
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In analogy to Eq. © we obtain the total relative sys- 
tematic error of P[n]/P[0]: 

(^)=|>„. r , + A„,(§M) (7) 

where A sys l denotes the systematic relative error of a 
single ratio. Both, the total and the individual systematic 
errors scale like the square of the statistical errors. Hence, 
if the statistical error is small the systematical error is 
negligible. By the same token, however, the method is 
not suitable for quickly generating a rough estimate of 
weight factors (to be used in a subsequent multicanonical 
simulation). If the statistical error becomes comparable 
to unity the systematic error will be of the same order. 

In Fig. we analyze the systematic error propagation 
for the runs with oj = 1 and 500,000 MC steps per state, 
which exhibit rather large statistical errors (cf. Fig. 3). 
Firstly, we regard the deviation between the probability 
distribution obtained from those short runs and an esti- 
mate for the true distribution obtained from a very long 
multicanonical run over the entire range. At small par- 
ticle numbers the systematic deviation agrees with the 
values calculated from Eq. © within the accuracy of the 
true estimate (0.5% — 2%). At large particle number, 
however, the calculated value is larger than the devia- 
tion from the true distribution. In this region, we sys- 
tematically underestimate the ratios ri upon increasing 
the number of particles, because we fail to equilibrate the 
configurations on the time scale of the short simulation 
runs. Therefore configurations with lower particle num- 
ber (contributing to Ha) are favored. These correlation 
effects are not included in Eq. Since both system- 
atical errors have opposite signs Eq. <0 constitutes an 
upper bound. 

In summary, we propose an extension to umbrella sam- 
pling in which simulation windows are sampled succes- 
sively. This allows us to generate weight factors by ex- 
trapolating results from previous windows. Hence, the 
scheme is able to calculate free energy barriers without 
prior knowledge of weight factors. The computational 
efficiency is competitive to alternative approaches which 
adjust weight factors on the fly and in principle indepen- 
dent from the choice of the window size. The algorithm 
is straightforward to implement and to parallelize. It in- 
volves neither fine-tuning of simulation parameters, nor 
violation of detailed balance. Errors can be calculated 
exactly and are controlled. A detailed analysis revealed 
the existence of a small systematic error which is, how- 
ever, irrelevant in practice. We also provided a scheme 
to check if configurations are well equilibrated. Our algo- 
rithm can be readily applied to a large variety of problems 
which involve the computation of free energy barriers. 



Possible applications include protein folding, the study 
of interfaces, nucleation barriers, and first order phase 
transitions. 
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